library(specr)
library(fixest)
library(readr)
library(broom.mixed)
library(texreg)
library(dotwhisker)
library(dplyr)
library(fixest)

# Two-Way linear Fixed Effect Models that use all data
# testing the effect of individual income sources 


source("funs_and_constants.R")

# add in data
epr_yearly_group_clusters_h1 <- read_csv("data/epro_data_h1.csv")

epr_yearly_group_clusters_h1$both_or_which <- relevel(factor(epr_yearly_group_clusters_h1$both_or_which), ref = "None")

## testing explicitly against each other
full_model_h1_cats <- feols(disagreement ~ both_or_which + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic
                        + n_tek_groups | year + countries_gwid, cluster = ~ countries_gwid,
                        data = epr_yearly_group_clusters_h1)

full_model_h1_cats_share <- feols(share_disagreement ~ both_or_which + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic
                                  + n_tek_groups | year + countries_gwid, cluster = ~ countries_gwid,
                                  data = epr_yearly_group_clusters_h1)

## cross-sectional
mean_h1 <- epr_yearly_group_clusters_h1 %>% 
  filter(year >= 2000) %>%# only obs after 2000 because of EPR 2.0 coding
  group_by(gwgroupid, countries_gwid) %>%
  summarize(
    across(c(share_disagreement, disagreement, n_orgs, regaut, groupsize, incidence_flag, status_pwrrank, any_multiethnic, nightlight_total_pc_log, n_tek_groups),
           ~ mean(.x, na.rm = T)),
    both_or_which = Mode(as.factor(both_or_which), na.rm = TRUE)
  )

full_model_h1_cats_cross <- feols(disagreement ~ both_or_which + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic + nightlight_total_pc_log + n_tek_groups | countries_gwid,
                       cluster = ~countries_gwid, data = mean_h1 %>% as.data.frame())

full_model_h1_max_cats_cross <- feols(share_disagreement ~ both_or_which + n_orgs + regaut + groupsize + incidence_flag +
                                        status_pwrrank + any_multiethnic + nightlight_total_pc_log + n_tek_groups | 
                                        countries_gwid,
                           cluster = ~countries_gwid, data = mean_h1)


screenreg(list(full_model_h1_cats_cross, full_model_h1_cats, full_model_h1_max_cats_cross, full_model_h1_cats_share),
          stars = stars,
          custom.coef.map = coef_map,
          custom.header = list("Disagreement (1/0)" = 1:2, "Prop. Disagreement" = 3:4),
          table = FALSE,
          include.proj.stats = FALSE)


texreg(list(full_model_h1_cats_cross, full_model_h1_cats, full_model_h1_max_cats_cross, full_model_h1_cats_share),
          stars = stars,
          custom.coef.map = coef_map,
          custom.header = list("Disagreement (1/0)" = 1:2, "Prop. Disagreement" = 3:4),
          table = FALSE,
          include.proj.stats = FALSE,
          file = "tables/h1_fourfold_comparison_robustness.tex")





## testing for 'direct' single-source effects, only nat resources is signficiant (results not in paper)
direct_nat_resources <- feols(disagreement ~ natural_resources_bin + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic
                              + n_tek_groups | year + countries_gwid, cluster = ~ countries_gwid,
                              data = epr_yearly_group_clusters_h1)

direct_agri <- feols(disagreement ~ agri_production_bin + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic
                     + n_tek_groups | year + countries_gwid, cluster = ~ countries_gwid,
                     data = epr_yearly_group_clusters_h1)

direct_pasture <- feols(disagreement ~ pasture_land_bin + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic
                        + n_tek_groups | year + countries_gwid, cluster = ~ countries_gwid,
                        data = epr_yearly_group_clusters_h1)





